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ABSTRACT 


A  technique  is  described  for  computing  internal  wave  in¬ 
duced  magnetic  field  and  magnetic  field  gradient  spectra  as 
these  would  be  measured  by  airborne  magnetic  sensors.  The  in¬ 
tended  application  is  remote  sensing  of  internal  wave  activity 
in  the  upper  layers  of  a  deep  ocean  in  the  wavelength  range 
less  than  1  km.  A  computer  program  (INTWAVE)  is  presented 
which  is  capable  of  accepting  data  from  towed  thermistor  chain 
measurements.  Numerical  results  obtained  with  INTWAVE  are  pre¬ 
sented  for  the  measured  internal  wave  spectra  reported  by  Nelson 
and  Milder  [M], 
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I.  INTRODUCTION 


The  subject  of  this  paper  is  the  determination  of  Internal 
wave  induced  magnetic  field  and  magnetic  field  gradient  spectra 
as  would  be  measured  by  a  magnetic  sensor  mounted  on  an  airborne 
measurement  platform.  The  principal  interest  is  in  sensing  of 
Internal  wave  activity  in  the  upper  layer  of  a  deep  ocean  in 
the  wavelength  range  below  1  km.  The  theoretical  basis  of  the 
approach  used  here  has  been  set  forth  in  [1]  in  conjunction  with 
the  internal  wave  model  developed  in  [3].  The  purpose  of  the 
present  paper  is  to  describe  a  systematic  technique  for  computa¬ 
tion  of  magnetic  field  and  magnetic  field  gradient  spectra  using; 
as  input  data  measured  "towed”  spectra  of  internal  wave  vertical 
displacement.  The  internal  wave  model  underlying  the  magnetic 
field  spectra  calculations  is  based  on  the  equipartition  hypoth¬ 
esis  of  Milder  [3], 

In  Sections  II  and  III  general  characteristics  of  Internal 
wave  spectra  predicted  by  this  model  are  summarized  and  the 
physical  constants  and  ocean  stratifications  considered  in  the 
calculations  are  presented.  Section  IV  provides  the  analytical 
relationships  between  the  internal  wave  parameters  and  the  mag¬ 
netic  field  and  magnetic  field  gradient  spectra  for  airborne 
observation  platforms. 

A  general  computer  program  (INTWAVE)  capable  of  providing 
the  spectra  for  rather  general  deep  ocean  Valsala  frequency  pro¬ 
files  is  described  in  detail  in  the  Appendix.  The  input  data 
can  be  either  in  the  form  of  analytic  expressions  or  supplied  in 
numerical  form. 


Representative  results  obtained  with  INTWAVE  are  plotted 
in  Pigs.  3-8.  Figures  3  and  4  give  magnetic  field  gradient 
spectra  for  the  exponential  thermocline  described  by  Garrett 
and  Munk  [2] .  Figures  5-8  are  based  on  the  measured  data  as 
reported  by  Nelson  and  Milder  [4],  the  VSisala  frequency  pro¬ 
file  for  which  is  plotted  in  Fig.  1  (labeled  case  4,  STD  data). 
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FIGURE  1.  Vaisala  Frequency  Profiles  Used  in  the  Calculation  of 
Magnetic  Field  Gradient  Spectra 
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II.  THE  INTERNAL  WAVE  MODEL 

The  model  adopted  herein  for  internal  waves  in  the  upper 
layers  of  a  deep  ocean  is  based  on  linearized  dynamics  under 
conditions  of  vanishing  mean  shear  flow.  The  statistical  de¬ 
scription  leading  to  analytical  representations  of  temporal 
and  spatial  spectra  utilizes  the  usual  assumptions  of  spatial 
and  temporal  statlonarity  as  employed,  e.g.,  by  Garrett  and 
Munk  [2],  The  interest  here  is  in  detection  of  internal  wave 
phenomena  in  the  frequency  range  u>  >  0.5  Nma^  (w  =  radian  fre¬ 
quency,  N  =  max  Vaisala  frequency),  ‘i.e.,  in  the  so-called 
max. 

buoyancy  range.  In  this  range  the  results  of  [2]  are  not 
directly  applicable. 

A  theory  which  focuses  on  the  buoyancy  range  has  been  ad¬ 
vanced  by  Milder  [ 3 J •  Milder 's  theory  is  based  on  rather  gen¬ 
eral  premises  Involving  conditions  of  statistical  equilibrium. 
The  underlying  hypothesis  is  that  of  equipartition  involving 
vorticity  and  strain  of  the  internal  wave  field.  The  conse¬ 
quences  of  this  equipartition  can  be  reduced  to  one  remarkably 
simple  rule:  the  energy  content  of  the  normal  modes  Is  propor¬ 
tional  to  the  square  of  the  mode  phase  velocities.  It  has  been 

shown  that  in  the  frequency  range  w  «  6 . 5  N  Milder 's  theory 

nicix 

gives  essentially  the  same  results  as  that  of  Garrett  and  Munk 
[3],  [1,  p.  265].  The  functional  form  of  the  spatial  spectrum 
of  the  isotherm  displacement  predicted  by  Milder • s  theory  ap¬ 
pears  to  be  in  agreement  with  towed  thermistor  measurement  data 
reported  by  Nelson  and  Milder  [4]. 

Specifically,  the  equipartition  hypothesis  leads  to  the 
following  canonical  form  of  the  spatial  wave  number  spectrum 
for  the  vertical  displacement: 


(1) 


KVK-y)  ‘  TFk-  (1  -  e'2Ky> 


where  y  (y  >  0)  Is  the  depth  and  K  =  2ir/A  is  the  horizontal  wave- 
number.  The  spatial  wavenumber  spectrum  S  (K,y)  Is  defined  such 
that  the  statistical  average  of  the  square  of  the  internal  wave 
vertical  displacement  is  given  by 


<  n  >  -  2n  J  dK  KSrin(K,y) 
0 


The  spectral  quantity  I(K)  depends  only  on  the  magnitude  of  the 
horizontal  wave  number  so  that  the  spatial  spectrum  S  (K)  is 
isotropic.  The  most  striking  feature  of  (1)  is  its  total  inde¬ 
pendence  of  the  ocean  stratification.  On  the  other  hand,  the 
spatial  spectra  of  the  fluid  velocity  components  depend  explic¬ 
itly  on  the  Vaisala  frequency  profile  N(y).  From  [1,  p.  255] 
the  spatial  spectrum  of  the  vertical  velocity  S  is  given  by 

y  y 

oo 

KSyy(K,y)  -  ^  I(K)  J  dy'  N2(y')  [e“K(y+y'  )  -  e“Kly“y'  G  2  .  (3) 


Similarly,  one  can  show  that  the  spectra  of  the  two  orthogonal 
horizontal  velocity  components  are 


00 

x(K,y)  -  — ^  I(K)  j'  dy '  N2(y')  [e_K(y+y,)  -  e"Kly-y'  0  , 


00 

KSzz(K,y)  -  S^I(K)  J"  dy'  N2(y» )  j  [e_K(y+y' 5  -  e_Kly"y' ']|2  .  (5) 
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The  spectra  of  the  components  of  the  horizontal  velocity  are 
not  isotropic  but  exhibit  the  characteristic  cos2w  (sin2w) 
dependence  on  w,  w  being  the  polar  angle  of  K.  On  the  basis 
of  data  in  [4]  one  finds  that  within  a  wide  •  nge  of  wavelengths, 
vi<=: . , 


60  m  <  X  <  1500  m  (6) 

the  spectral  excitation  function  I(K)  has  the  form 

I(K)  =  Cp  K“p  ,  (7) 

where  Cp  is  a  constant  dependent  on  the  exponent  p.  The  defini¬ 
tions  of  I ( K )  in  this  paper  and  [1]  differ  from  the  definition 
in  [3].  If  we  denote  the  I(K)  given  in  [3,  Eq.  52,  p.  22]  by 
I ( K ) ,  then 


I (K)  =  i  I(K)  .  (8) 

8iTd 

Also,  the  spectral  density  S._  of  the  vertical  displacement  as 
given  by  Eq .  (52)  in  [3]  is  .related  to  the  spectral  density 
in  (1)  by 


S 


on 


(2tt) 


(9) 


The  relationship  (8)  can  be  inferred  by  comparison  of  Eq.  (5) 
on  page  5  of  [3]  with  Eqs.  E-71  and  E-73  of  [1], 

In  the  calculations  presented  in  Section  IV  we  consider 
three  types  of  stratifications.  These  are  shown  in  Pig.  1. 

The  sharply  peaked  profile  identified  in  the  figure  as  "case  4  - 
STD  data  set"  is  taken  from  Pig.  3  of  [3],  which  is  based  on 
the  towed  thermistor  chain  measurements  of  internal  waves  in 
the  eastern  tropical  Pacific  Ocean  in  the  early  fall  of  1974, 
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as  reported  in  [4].  The  real  data  extends  only  to  about  220 
meters.  For  the  purposes  of  analytical  modeling,  an  exponen¬ 
tially  decreasing  VaisalS  frequency  profile  of  the  form 
exp  (-  y/1300)  (y  in  meters)  is  assumed  for  y  >  220  meters. 

The  profiles  labeled  case  1  and  case  2  are  those  used  by  Garrett 
and  flunk  [2].  In  particular,  in  the  exponentially  stratified 
ocean  model  (case  1),  N(y)  -  5*236  x  1CT3  exp  (-  y/1300)  rad/sec, 
i.e.,  no  mixed  layer  is  assumed.  In  case  2,  a  mixed  layer  is 
assumed  to  a  depth  of  50  m  (i.e.,  N(y)  “0,  0  s  y  <  50  m) ;  for 
y  >  50  m  the  same  exponential  decay  law  is  assumed  as  in  case  1. 

For  the  STD  data  set,  the  empirically  determined  constants 
C  and  p  are  given  in  Eq.  (5*0  on  page  23  of  [3].  The  units 
employed  in  [1]  and  in  the  magnetic  field  calculations  in  Sec¬ 
tion  IV  are  always  m,  kg,  sec.  In  particular,  the  excitation 
function  I(K),  as  given  by  (1)  (and  also  in  [1])  always  has 
dimensions  of  m  ,  while  the  wavenumber  K  =  2tt/A  .  In  Eq.  (5*0 
of  [3]  different  units  are  employed.  In  addition,  private 
communication  with  the  author  of  [3]  revealed  that  k  in  Eq.  (5*0 
of  [3]  is  actually  1/A  instead  of  2-n/X  (A  is  in  km).  Taking 
this  into  account,  as  well  as  the  different  normalization  of 
the  excitation  function  as  shown  in  (8),  the  constant  Cp  to  be 
used  in  (7)  is 


C 


1 


5 


Consequently, 


.12  x  (2ti)2'5 
101*5  x  8tt2 


4.756  x  10"3  (meters)1/2  . 

(10) 


I (K)  =  4.756  x  IQ-3  K~1,5  m2  ■  .  (11) 

For  the  exponential  profiles,  Cp  is  determined  from  the 
energy  density  constant  E  and  cutoff  wavenumber  given  in  [2]. 

The  cutoff  wavenumber,  Kc ,  corresponds  to  the  lower  truncation 
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limit  in  the  Garrett  and  Munk  wavenumber  spectrum.  We  express 
it  in  terms  of  the  dimensionless  quantity  v  , 

V 

vc  =  Kflb  =  .3267  (12) 

where  b  is  the  depth  at  which  the  Vaisala  frequency  has  decayed 
to  1/e  of  its  maximum  value  (at  the  surface  or  just  below  the 
mixed  layer).  The  second  fundamental  constant  is  E,  the  inte¬ 
grated  energy  density  c  .-r  a  water  column.  The  value  given  in 
[2]  is 

E  -  3-82  x  103  joules/m2  .  (13) 


The  two  exponential  profiles  in  Fig.  1  can  be  represented  ana¬ 
lytically  as  follows: 


(°,  0  <  y  <  yB 

H  exp  {(3vy)/b}  ;  y  >  om 


where  =  0  and  ym  =  50  m  correspond,  respectively,  to  case  1 

and  case  2.  The  constant  Cp  is  evaluated  by  integrating  the 

spatial  energy  spectrum  between  K  =  K  and  K  B  »,  as  described 

c 

on  page  266  of  [1].  For  the  general  case  in  (14)  the  formula 
for  C  becomes 

Jr 


no  N^bp+1f  (y  /b) 
Ho  m  pvjrnr  ' 


(15) 


where 


/*“  r  -2y_v/b  l 


f  (y  /b ) 

p  'm 


(16) 
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3  2 

where  pQ  =  10  kg/m  and  r  Is  the  Gamma  function.  For  y^  ■  0, 

*  one  finds  [1] 

C2  -  1.2166  x  10"3 

l  *  C2  -  1.8727  x  10-2  m 

f  When  ym  /  0  the  result  of  the  integration  in  (16)  cannot  be  ex¬ 

pressed  in  terms  of  tabulated  functions.  A  numerical  integra- 
:  *  tion  subroutine  is  incorporated  in  the  general  computer  program 

in  the  Appendix. 

Garrett  and  Munk  use  p  B  2  in  their  paper  [2],  In  order  to 
show  the  sensitivity  of  the  results  to  changes  the  two  extreme 
values  p  ■  1  and  p  =  2  will  be  considered. 
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III.  RELATIONSHIP  BETWEEN  EXPERIMENTALLY  DETERMINED 
TOWED  SPECTRA  AND  THE  EXCITATION  FUNCTION  I ( K) 


» 


I 
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References  [3]  and  [4]  provide  plots  of  spatial  spectra 
at  different  depths  from  which  the  constants  given  in  (11)  are 
determined.  These  spectra  are  identified  as  "towed  spectra", 
presumably  obtained  from  an  analysis  of  time  series  as  measured 
by  thermistors  towed  through  the  water.  The  detailed  steps  em¬ 
ployed  in  the  process  of  reducing  the  measured  time  series  to 
the  spatial  spectra  are  not  presented  in  [3]  and  [4].  There 
are  several  approximations  associated  with  the  translation  of 
towed  time  series  data  to  equivalent  spatial  spectra,  which 
approximations  also  have  bearing  on  the  interpretation  of  the 
results  for  magnetic  field  spectra  presented  in  Section  IV. 

The  following  Is  a  brief  discussion  of  the  nature  of  these 
approximations . 

If  we  denote  the  Instantaneous  vertical  displacement  at 
depth  y  in  the  water  by  n(y,£,t)  (where  p  denotes  the  horizon¬ 
tal  coordinates),  the  temporal  spectrum  0^^(y,u)  observed  in 
the  water  at  a  point  moving  with  horizontal  velocity  U  may  be 
defined  formally  by  the  relationship 


«'u)  <y.»>  - 
nn 


lim 

T 


< 


/ 

-T/2 


T/2 


e~lu)t  n(y,Ut,t)  dt 


> 


(17) 


In  an  actual  measurement  T  would,  of  course,  be  finite,  as  a  re¬ 
sult  of  which  the  spectral  resolution  would  be  limited  to  1/T. 
The  temporal  spectrum  4>^^(y,u>)  cannot,  in  general,  be  expressed 
in  terms  of  the  spatial  spectrum  S  (K) .  This  can  be  done, 
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however,  approximately  for  fast  tow  speeds,  In  which  case  the 
following  asymptotic  result  is  obtained  [1]: 

it/2 

^Ky)  ~  jj  J  d<3  I(a)U-1sec0)  (1  -  exp  -  2yoU”1sece)  ,  (18) 

0 

where  I  is  the  excitation  function  entering  in  the  definition 
of  the  spatial  spectrum  in  (1).  Equation  (18)  holds  for  U  much 
greater  than  the  group  velocity  of  the  internal  wave  motion 
(fast  tow  speed  approximation).  In  addition,  (18)  is  not  valid 
for  w  at  and  below  the  maximum  VaisalM  frequency.  When  (18)  is 
evaluated  at  to  ■  KU,  one  obtains 

tt/2 

^nn5  (KU,y)  ~  O’  f  d(3  --p--  (1  -  exp-2Kysec0)  .  (19) 

0 

Apart  from  a  constant  factor,  (19)  is  identical  to  the  spatial 
spectrum  along  a  typical  "cut”  through  the  isotropic  internal 
wave  field.  In  [3]  and  [^]  the  form  of  I(K)  is  inferred  from 
plots  of  (KU,y )  as  a  function  of  K  at  various  depths  y. 

In  general,  the  determination  of  T(K)  would  require  solving 
the  integral  equation  in  terms  of  the  given  data.  If,  however, 

I ( K )  Is  aesurned  to  obey  a  power  law,  viz.,  (7),  then  (19)  becomes 
c  K~P-1  tt/2 

^^(KU,y) - -  C  cosp8  [1  -  exp-2Kysecg]  d3  (20) 

0 

For  sufficiently  large  depths,  viz.,  when 

2Ky  >>  1,  (21) 

the  exponential  term  in  the  integrand  may  be  dropped,  while 
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v/2  r  fv+1\ 

f  cosPfidB  -  ^  V 2“ ) 

J  cos  m  “  T  /p+2\ 


Consequently,  the  towed  spectrum  becomes 


U4^}(KU,y) 


r(P+l) 

t 'l \  2  /  r  .,-p-l 

2  „  /p+2\  CpK 


permitting  the  empirical  determination  of  C  and  p  from 
^(KU.y).  P 


(22) 
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IV.  MAGNETIC  FIELD  COMPONENT  AND  GRADIENT  SPECTRA 
OBSERVED  FROM  AIRBORNE  MEASUREMENT  PLATFORMS 

We  assume  a  platform  moving  with  velocity  V  m/sec  at  a 
fixed  height  of  y  above  the  ocean  surface  (Pig.  2).  The  geo¬ 
magnetic  field  is  in  the  xy  plane.  Its  direction  with  respect 
to  the  horizontal  plane  (sea  surface)  is  defined  by  the  angle 
<J>D  (dip  angle).  The  motion  of  the  platform  is  parallel  to 
which  is  oriented  at  an  angle  a  with  respect  to  the  vertical 
plane  containing  the  geomagnetic  field.  The  components  of  the 
magnetic  field  induced  by  internal  waves  are  resolved  along  the 
mutually  perpendicular  directions  4^,  which  may  be  con¬ 

sidered  fixed  with  respect  to  the  platform.  The  temporal  cross¬ 
spectrum  of  magnetic  field  components  Bv,  will  be  denoted  by 
S  .  A  typical  component  of  the  magnetic  field  gradient,  G^, 
is  defined  as  the  derivative  of  B  with  respect  to  direction  v. 
The  gradient  components  are  symmetrical,  i.e.,  GV(j  =  Gyv' 
temporal  cross  spectra  of  the  gradients  are  denoted  by  SVy . . 

V 


1)  OCEAN  SURFACE  PARALLEL 

TO  U  PURE 

2)  DIP  ARSLE 

I)  OEOMMRETIC  FIELD  IR 
xy  PURE 

4)  INDUCED  MA8RETIC  FIELDS 
ARD  6RADIERTS  RESOLVED 
AL0R6  URIT  VECTORS 

-I'^l'-S 

8)  PLATFORM  MOTION 
AL0H0i.i 
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For  frequencies  above  the  maximum  Vaisala  frequency,  the 
magnetic  field  and  gradient  spectra  can  be  approximated  by  the 
following  asymptotic  expressions  [1]  (large  tow  speed  approxi¬ 
mation)  : 


(oy  B)2  C.  71/2 

Svy(u),y)~  ““  “72 


2—  f  dp  g^  (a,3)  K"p+1e"2Ky  /'y'2N2(y')e2Ky'dy'  , 

2(1+  3  cos"  <j)D>  £  ^  £ 


(23) 


Svyivy(w»y)< 


Cp(c™o  B) 


2  tt/2 


(1  +  3  cos^  <j>D)u>  £ 


/  d6  ®vy;vy(oi,0)  K ~P+3e“2KyJ' y'2N2(y' )e2Ky’dy’ 

—00 

(24) 


where  y  is  the  height  above  the  ocean  in  meters,  K  »  w/(Vcos3), 
o  =  4  mho/m,  y  =  4tt  x  10“^  h/m,  B  ■  6.24  x  10^  pT  and  the  con- 

0  O 

stant  C  has  dimensions  of  (meters)  p.  The  quantities  g  in 
the  integrands  are  dimensionless  trigonometric  polynomials, 
defined  on  page  A-l4.  The  spectra  in  (23)  and  (24)  have  dimen- 
sions  of  (pT)  /Hz  and  (pT/m)  /Hz,  respectively.  The  Internal 
wave  vertical  displacement  spectrum  associated  with  the  magne¬ 
tic  field  spectra  in  (23)  and  (24)  is  given  by  (20).  As  noted 
earlier,  according  to  the  equipartit ion  hypothesis,  the  latter 
is  Independent  of  N(y).  On  the  other  hand,  the  magnetic  fields 
are  seen  to  depend  explicitly  on  N(y),  (C.f.  Eqs .  3-5.) 

The  computer  program  described  in  the  Appendix  is  capable 
of  providing  magnetic  field  and  gradient  spectra  for  an  arbi¬ 
trary  set  of  parameters  and  Vaisala  frequency  profiles.  The 
N(y)  profile  input  data  can  be  supplied  either  analytically  or 
as  a  set  of  values  of  N(y)  spanning  a  specified  depth  below  the 
ocean  surface.  Below  this  depth  an  exponential  decay  law  is 
assumed.  The  constants  and  p  are  to  be  supplied  with  the 
data. 
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Representative  results  obtained  with  the  aid  of  the  INTWAVE 
computer  program  are  plotted  in  Figs.  3-8.  The  direction  of 
platform  motion  a  (defined  in  Fig.  2),  platform  velocity  V 
(m/sec),  height  above  the  ocean  surface  Y  (rn),  the  magnetic 
dip  angle,  and  the  exponent  p  governing  the  decay  of  the 
wave  number  spectrum  (Eq.  7)  are  shown  in  each  figure.  The 
spectra  are  in  (pT/m)  /Hz,  the  horizontal  scale  Is  in  Hz,  and 
the  subscripts  pertaining  to  the  gradient  components  conform 
to  the  coordinate  system  Z^,  >  i-n  Fig.  2. 

Figure  3  shows  the  vertical-horizontal  (along  the  track) 
gradient  spectrum  for  the  exponentially  stratified  ocean  with¬ 
out  a  mixed  layer.  Figure  4  shows  the  spectrum  for  the  same 
gradient  component  for  the  exponentially  stratified  ocean  with 
a  mixed  layer  (case  2,  in  Fig.  1).  The  differences  between 
Fig.  3  and  Fig.  4  are  seen  to  be  quite  insignificant.  Figure  5 
shows  the  spectra  of  the  three  orthogonal  gradient  components 
for  the  STD  data  set.  Although  the  functional  form  of  the 
spectra  is  essentially  the  same  as  for  the  exponentially  strat¬ 
ified  ocean,  the  overall  levels  of  the  magnetic  field  gradients 
are  substantially  higher.  This  is,  of  course,  to  be  expected, 
in  view  of  the  much  higher  peak  Vaisala  frequency  in  this  case. 
Somewhat  lower  spectral  values  are  obtained  when  the  platform 
motion  is  oriented  at  90  deg  with  respect  to  the  plane  of  the 
geomagnetic  field,  as  shown  in  Fig.  6.  In  Figs.  3-6  the  dip 
angle  is  taken  as  zero  (equatorial  zone).  A  dramatic  Increase 
in  level  is  obtained  when  the  dip  angle  is  taken  as  90  deg 
(polax'*  regions),  as  shown  in  Fig.  7  for  the  STD  data  set. 

I 

In  Figs.  3-7  the  platform  velocity  is  100  m/sec  (typical 
aircraft  speed).  To  Illustrate  the  effect  of  platform  velocity, 
Fig.  8  shows  the  vertical-horizontal  (along-the-track)  gradient 
spectrum  for  platform  speeds  of  100  m/sec  and  10  m/sec. 
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Gradient  Spectra 
stratified  ocean) 


(pT/m)2/Hz 


FIGURE  4.  Internal  Wave  Induced  Magnetic  Field  Gradient  Spectra 
from  a  Moving  Platform  (exponentially  stratified  ocean 
with  a  mixed  layer  y  ■  50  meters) 
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Internal  Wave  Induced  Magnetic  Field  Gradient  Spect 
from  a  Moving  Platform  (STD  Data  Set) 


(pT/m)2/Hz 


FIGURE  6.  Internal  Wave  Induced  Magnetic  Field  Gradient  Spectra 
from  a  Moving  Platform  (STD  Data  Set) 


(pT/Bi)2/Hz 


FIGURE  8.  Internal  Wave  Induced  Magnetic  Field  Gradient  Spectra 
from  a  Moving  Platform  (STD  Data  Set) 


V.  CONCLUSIONS 


A  general  conclusion  that  can  be  drawn  from  the  results 
presented  herein,  as  well  as  from  the  results  of  a  large  number 
of  similar  calculations  using  the  INTWAVE  program,  is  that  in¬ 
ternal  wave  induced  magnetic  field  gradient  spectra  observed 
from  a  moving  aircraft  display  no  interesting  structural  fea¬ 
tures:  all  the  spectra  decrease  monotonically  with  increasing 
frequency  and  are  quite  insensitive  (both  in  regard  to  level 
and  functional  form)  to  the  spatial  variation  in  the  Vaisala 
frequency  profile.  The  critical  parameters  that  affect  the 
spectral  level  appear  to  be  the  area  enclosed  by  the  N(y) 
profile  and  the  spatial  wavenumber  decay  constant  p,  as  well 
as  the  geomagnetic  latitude. 
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I.  PURPOSE 

The  program  INTWAVE  was  designed  to  calculate  and  plot 
the  Internal  wave-induced  magnetic  field  spectra  (component 
and/or  gradient)  from  a  moving  platform  as  a  function  of  fre¬ 
quency.  It  provides  for  a  choice  among  various  assumptions 

■  »  H  •• 

as  to  the  behaviour  of  the  VAISALA  frequency  vs.  ocean  depth. 

II.  EQUATIONS 


The  expression  for  the  magnetic  field  component  Bpeotrum 
is  given  by  Equation  255  of  Reference  1.  Rewritten  in  the 
notation  of  the  INTWAVE  program,  the  equation  becomes: 

tt/2 

-p+1  2my 


T1  I  «*<'»>  (wgis)" 


S±i  a  C0^ST1  /  G01(e)  [ I  e  Vcosp  FN(S)  d0  (1) 


where 


C0NST1 


49190 


l  +  ScOS^Q 


(2) 


f  2u>jf 

PN(e)  “  Cp  I  Y2N2(Y)  eVcoB$  dY, 


and 


(3) 


aci(e)  -  G^Ca+B-ir)  +  (oc-3-tr )  (4) 

The  expressions  (Equation  4)  are  a  set  of  trigonometric 
polynomials  defined  for  i,j  ■  1,3.  These  polynomials  are  given 
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in  Addendum  1.  The  N(Y)  in  Equation  3  is  the  VAISALA  .frequency 
for  ocean  depth  Y ,  and  Cp  is  an  associated  constant  supplied 
with  each  N(Y)  function. 

Similarly,  the  expression  for  the  magnetic  field  gradient 
spectrum,  given  by  Equation  259  of  Reference  1  becomes  in  the 
program  notation 


kJlmn 


2C0NST1 

K  »  ■■  ■■"■.■■i..  — 

(D 


— ion  - 

e  Vcos$  pN( 3 )  dp 


(5) 


where 


C0NST1  and  PN ( B )  are  defined  as  for  Equation  1  and 

0C2<»>  ■  +  <W 

The  Q^^mn  expressions  of  Equation  6  are  a  Bet  of  trigono¬ 
metric  polynomials  defined  for  i,J,k  *  1,3  and  are  listed  in 
Addendum  1. 


III.  PARAMETERS 


The  parameters  which  are  required  for  a  particular  cal¬ 
culation  are  as  follows:* 

1.  V  -  platform  veloolty  (m/sec) 

2.  y  -  height  of  platform  above  ocean  surface  (m) 

3.  a  (Alpha)  -  angle  between  the  direction  of  platform 
motion  and  the  plane  of  the  geomagnetic  field  (rad.) 
(0<ci<_27t) 

4.  4>[)TPhiD)  -  magnetic  dip  angle  (rad.)  ( 0</|)D<_7r/2 ) 

5.  i(ii)  and  j(v1J)  -  directions  of  components  of  the  in¬ 
duced  magnetic  field  (i,j*l,3) 

If - 

The  program  name  for  each  variable,  when  different,  is  given 
in  parenthesis.  In  particular,  all  Greek  letters  are  spelled 
out  in  the  program. 
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p»wni'»i'ii 1  . -  Trv-y.  - 


6.  k(kk),  &(HJl),  m(mm),  n(nn)  -  the  directions  of  gradi¬ 
ents  of  the  induced  magnetic  field  (k,Jl,m,n  -  1,3;  k<&) 

7.  Ncase  -  an  integer  associated  with  each  assumption 
about  the  VAISALA  frequency  vs.  ocean  depth  (see  Sec¬ 
tion  IV)  (l<Ncase<4) 

8.  p  -  the  internal  wave  mode  decay  constant.  For  Ncase* 
1,2,  or  3,  P  ■  1  or  2.  For  Ncase«4,  p>0. 

9.  Yra(YM)  -  this  parameter  is  required  only  when  Ncase* 

2  or  3.  It  is  the  value  of  the  depth  at  which  the  peak 
value  of  N(Y)  occurs  and  is  restricted  to  0<Ym<200m 

The  calculations  of  S^j  and/or  Sk^mn  are  made  as  a  func¬ 
tion  of  w  (omega),  over  a  range  from  in  *  .02ir  to  2ir,  the  values 
being  more  closely  spaced  at  the  low  end  and  more  widely  spaced 
near  2u>.  The  actual  spacing  is  given  by  u/2w  *  .01(.0l).l; 

. 1 ( . 02 ) . 2  ;  ,2( .025)1.00. 

IV.  THE  VAISALA  FREQUENCY  VS.  DEPTH  OF  OCEAN 


The  function  FN(ft)  in  liquation  3  is  dependent  on  the  na¬ 
ture  of  the  function  N(Y),  where  Y  Is  ocean  depth  (m),  (Y<0) 

M  M  •» 

and  N ( Y )  is  the  VAISALA  frequency  for  that  depth.  Several  pos¬ 
sible  assumptions  as  to  the  nature  of  this  relationship  have 
been  Incorporated  into  the  program,  Including  one  which  is 
general  In  nature  and  can  accept  user  supplied  data  points. 

They  all  contain  a  parameter  b  which  is  the  reciprocal  of  the 

•  •  •  ■  it 

VAISALA  frequency  asymptotic  decay  constant.  In  the  program, 
at  present  the  value  of  b  is  set  at  1300  meters,  but  could 
easily  be  changed. 

The  four  different  assumptions  which  the  program  Is  de¬ 
signed  to  handle,  indicated  by  the  value  of  the  input  parameter 
Ncase,  are  as, follows: 

1.  Ncase  =  1  assumes  an  exponentially  stratified  ocean 
given  by 

N (Y )  -  Nm  eY/b  (7) 

When  this  Is  substituted  in  Equation  3  and  the  re 
suiting  expressions  Integrated,  we  obtain 
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MMHWH  ii  1  'Jrr.pa*  ■ 


where 


PN($)  =  Kp 


oobB 


|oob8+  tub  I 


(8) 


K±  =  .2169b 

K2  -  .1832 


2.  Ncasei=2  assumes  an  exponentially  stratified  ocean 
with  a  mixed  layer.  N(Y)  is  defined  by: 


N(Y)-0  ,-Ym<  Y<0 

N(Y)-Nme(y+Ym^/b  ,  Y<~Ym 


(9) 


The  value  of  Ym  Is  user  supplied. 

Substituting  this  expression  In  Equation  3»  and  inte¬ 
grating  the  result  yields 


HU)  -  Kp.eVcos3-f(B)* 


f ( 3)  +f2(B) 


(10) 


where 


f  (S) 


cosB 

oosB+  mb 
V 


V  „  .  3040b 

Ki  '  WmJ^TD 


V  m  .3040 

2  vTnT2TY^2T 


and 


VINT2(Yra,p) 


-2YmV/b 

1+V 


dV 


(11) 


where 


Vc  -  .3267 

The  values  of  VINT2(Yra,p)  were  computed  by  Gaussian 
numerical  quadrature  for  p"l,'1  and  for  ym  “  0(25)200.  The 
results  were  such  that  linear  interpolation  for  Ym  is  adequate. 


3.  Nca3e=3  assumes  a  thermocline  with  a  well-defined  peak 
at  Ym.  N(Y)  is  given  by: 

N(Y)  =  Nm(?|)Ym/b  e(Y+Ym)/b  (12) 

Integrating  the  resulting  expression  for  PN(0)  yields 

,  .  (- T ^  +  3 ) 

PN(S)  -  Kp  +  l)  jcooeL)  (13) 

where 

K1  “  V I  NT  3  TH  and  K2  “  VINT  3^  ,27 

CO 

VXNT3  ( Ym,p )  -  J  V-P-1  |l  -  av  (I'D 

Vc 

The  values  of  VINT3(Ym,p)  were  computed  in  the  same  way 
and  for  the  same  values  of  Ym»p  as  for  VINT2(Ym,p).  Again 
linear  Interpolation  for  0<Ym<200  is  adequate. 


NOTE:  If  Ym“0,  the  previous  three  cases  become  identical. 

For  cases  2  and  3,  a  user  supplied  value  of  Ym<200  should  be 
supplied. 

4.  Ncase=4.  This  is  the  most  general  case  in  that  the 
N(Y)  vs.  depth  relationship  up  to  some  maximum  depth 
of  YMIN  is  described  by  a  user  supplied  set  of  up  to 
100  data  points. #  For  depths  below  YMIN,  an  exponen¬ 
tial  relationship  is  assumed. 

A  set  of  NPT  pairs  of  points  NVAL ,  YVAL  from  depth  of  0 
to  depth  YMIN  is  supplied  where  NVAL«N(Y)/Nm,  so  that  0<NVAL<1. 


s - - — - — 

This  limit  on  the  number  of  data  points  could  easily  be 
changed  if  more  points  were  necessary. 


and  YVAL=depth,  so  that  0 _>YVAL>YMIN .  Since  any  value  of  N(Y) 
for  Y>YMIN  is  obtained  by  linear  interpolation  in  this  set  of 
points,  it  should  be  noted  that  the  points  should  (if  obtained 
from  a  curve)  be  chosen  in  such  a  way  that  linear  interpolation 
is  satisfactory,  i.e.,  points  should  be  more  closely  spaced 
where  curvature  is  greatest  and  should  include  all  minima  and 
maxima.  If  YVAL  =>  NMIN  at  Y  *  YMIN,  the  exponential  expression 
used  for  N(Y)  below  depth  YMIN  is 


N(Y)  -  N  NMIN  e(y“YMIN)/b 
'  '  m 

The  expression  for  FN(g)  for  Ncase*4  becomes 

/  } 

FN(3)  =  Cp  l  J  Y2N2(Y)eVcosB  dY 
(YMIN 

2aYMIN 

+  Nm?-NMIN2eVcOBE  [(aYMIN-l)2  +  l] 
a^  J 


(15) 


(16) 


where 


a 


2(i 


_j<i  \ 
Vcbs 


The  Integral  from  Y  “  YMIN  to  0  Is  obtained  by  numerical 
integration  (Gaussian  quadrature  formula). 

In  using  Ncase=i(,  along  with  the  set  of  data  points  YVAL, 
NVAL,  the  set  of  parameters  NBet  (an  identifier),  NPT,  Nm, 
NMIN,  YMIN,  and  Cp  must  all  be  supplied  by  the  user.  For 
Nset=l,  a  set  of  STD  data  points  and  parameters  has  been  used, 
the  values  of  which  are  given  in  Addendum  2. 
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V.  NUMERICAL  QUADRATURE 


A.  THE  sjj»  SkHnin  ORATION 

The  integration  indicated  by  Equations  1  and  5  is  per¬ 
formed  in  this  program  by  using  SIMPSON's  rule  as  follows. 

The  integral  is  obtained  by  using  the  statement 

INT  =  SIMP 2 ( BEGIN, END, ACC ,NN,FUNCT,NN1 , ACC 1 ) 
where  the  input  parameters  are 

FUNCTN  -  The  integrand,  defined  in  a  function 

subroutine 

BEGIN  -  The  beginning  of  the  interval  of  inte¬ 

gration  (  =  0) 

END  -  The  end  of  the  interval  of  Integration 

(  =  ir/2 ) 

ACC  -  The  accuracy  desired  in  the  result 

(=.0001) 

NN  -  The  maximum  number  of  intevals  into 

which  the  interval  is  to  be  subdivided 
before  terminating  the  iteration  (-1000) 

and  the  outputs  are 

NN1  -  Actual  number  of  intervals  before  itera¬ 

tion  was  terminated 

ACC1  -  Relative  accuracy*  achieved  in  final 

estimate 

INT  -  The  final  estimate  of  the  value  of  the 

integral 


If  the  final  estimate  is  AREA2  and  the  previous  estimate  is 
AREA1,  the  relative  accuracy  Is  defined  as  AREA2-AREA1 


An  initial  estimate  of  the  integral  is  obtained  by  using 
four  intervals  and  then  the  number  of  intervals  is  successively 
doubled  until  either  NN1>NN  or  ACC1<ACC. 

In  actual  use  for  many  combinations  of  parameter  inputs, 
it  was  found  that  in  most  cases  only  16  or  32  intervals  were 
required  to  achieve  accuracy  of  .0001. 

B.  INTEGRAL  IN  FN ( B)  for  Case  4 

The  integral  in  Equation  16  must  be  obtained  by  numerical 
quadrature  since  N(Y)  vs.  Y  is  given  by  a  tabular  set  of  values, 
linear  interpolation  being  used  to  obtain  intermediary  values. 

An  attempt  to  get  accuracy  of  .0001  with  SIMPSON’s  rule  re¬ 
quired  on  the  order  of  2000-8000  points  per  integration;  there¬ 
fore,  the  Gauss  method  of  numerical  Integration  (in  which  the 
points  of  subdivision  are  not  equally  spaced,  but  chosen  to 
give  greater  accuracy)  was  substituted  using  96  points. 

VI.  INPUT  CARD  DECK 

This  section  provides  the  information  needed  to  supply  an 
input  deck  required  for  a  computer  run  of  INTWAVE. 

With  the  exception  of  the  YVAL,  NVAL  set  of  data  points 
supplied  when  Ncase-4,  all  inputs  are  members  of  a  NAMELIST 
group,  and  can  thus  be  given  values  by  means  of  a  NAMELIST  in¬ 
put  card  in  an  unformatted  manner,  the  details  of  which  are 
given  in  Addendum  3 • 

a.  The  three  NAMELISTS  which  contain  all  the  required 
Inputs  are: 

1.  OPTIONS  namelist.  This  is  a  group  of  five  option 
parameters  which  are  to  be  set  equal  to  0  or  1. 

The  non-zero  ones  must  be  defined  on  the  OPTIONS 
card.  The  parameters  are: 

(a)  IPRINT,  This  is  set-1  when  an  extra  printout 
is  desired  of  17  equally  spaced  values  of  the 
integrand  (from  B-  0  to  3  -  tt/2)  of  Equations 
1  and  5.  This  33  normally  used  only  in  a  debug 
run. 
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(b)  ICALCG.  This  is  set=l  if  values  of  the  component 
spectrum  (S-y)  are  to  be  calculated. 

(c)  ICALCC.  This  is  set=l  if  valueB  of  the  gradient 
spectrum  (S^^mn)  are  to  be  calculated. 


(d)  IPLOTC.  This  is  set=l  if  plots  of  the  component 
spectrum  (sij )  are  to  be  produced. 

(e)  IPLOTG.  This  is  set=l  if  plots  of  the  gradient 

spectrum  are  to  be  produced. 


WARNING  NOTES.  Plots  of  either  type  cannot  be  produced  unless 
the  corresponding  calculations  are  made,  whereas  the  inverse 
is  not  true.  Thus,  the  combination  ICALCG=0,  IPL0TG=1  Is  not 
a  feasible  one.  Also  when  IPRINT»1,  all  six  subscripts  ii  to 
nn  must  be  assigned  non-zero  values. 

2.  PARAMS  namelist.  This  group  contains  all  the  para¬ 
meters  listed  in  Section  III,  namely  V,  y,  a,  <J>d, 
p,  Ncase,  Ym,  and  the  six  subscripts  Ii  to  nn. 

(a)  Of  these  inputs  Ncase  and  the  six  subscripts 
are  integers;  the  others  should  have  a  decimal 
point  Included  in  their  value. 

(b)  Ym  is  needed  only  if  Ncase  =  2  or  3. 

(c)  The  subscripts  ii,jj  are  needed  only  if  ICALCC^l 
or  IPRINT-1. 

(d)  The  subscripts  kk,  Hi,  mm,  nn  are  needed  only 
if  ICALCG=1  or  IPRINT-1. 

(e)  Any  parameter  omitted  from  the  PARAM  card  will 
be  assigned  a  value  0,  unless  it  follows  a 
previous  PARAM  card  in  which  case  Its  value 
will  remain  unchanged. 

3.  CASE4  namelist.  This  set  of  inputs  is  required 
only  when  Ncase=4  on  the  PARAM  card  preceding  it. 

It  contains  values  for  the  following  variables, 
all  of  which  are  associated  with  the  set  of  N(Y) 
data  points  supplied. 

(a)  Nset.  This  is  an  integer  number  assigned  to 
each  set  of  data  points  to  uniquely  Identify  it. 

(b)  NPT.  The  number  of  data  points  supplied  for  the 
arrays  NVAL ,  YVAL . 

(c)  NMIN .  The  value  of  N(Y)/Nm  at  depth  YMIN. 

(d)  YMIN.  The  last  tabulated  depth  value  (a  nega¬ 
tive  number). 
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(e)  Nra.  The  maximum  value  of  N(Y). 

( f )  Cp  . 

NOTE :  Only  the  first  two  of  these  five  variables  are  integers. 

b.  CASE4  Data  Points.  The  CASE4  namelist  input  card 
must  be  immediately  followed  by  a  set  of  data  cards 
which  provide  a  set  of  NPT  values  for  the  NVAL,  YVAL 
arrays.  The  first  set  of  values  should  be  for  Y=0 
and  the  last  set  for  Y=YMIN.  The  array  of  positive 
numbers  -  YVAL  is  punched  on  the  input  cards  and  con¬ 
verted  in  the  program  to  negative  values.  The  format 
in  which  these  cards  are  to  be  punched  is  as  follows: 

1.  Each  pair  of  values  NVAL,  YVAL  occupies  a  field  of 
10  columns,  the  first  five  for  NVAL,  the  second 
five  for  YVAL. 

2.  The  number  may  be  placed  anywhere  within  the  5- 
column  field  as  long  as  the  decimal  point  is  in¬ 
cluded  . 

3.  Eight  pairs  of  values  are  included  on  each  card, 
until  all  values  are  punched. 

£.  Varying  Parameters  Within  One  Run.  A  complete  input 
deck  for  one  run  with  one  parameter  set  consists  of 
an  OPTIONS  card  followed  by  a  PARAMS  card,  followed, 
if  Ncase=4,  by  a  CASE4  card  and  associated  NVAL, YVAL 
cards,  all  terminated  by  an  end-of-file  card. 

Results  may  be  obtained  for  various  parameter  sets 
within  one  run  by  following  the  first  sequence  of  in¬ 
put  cards  by  a  string  of  PARAM  cards  before  the  end- 
of-file  card. 

Everytime  a  calculation  is  completed,  the  next  PARAM 
card  in  the  input  deck  is  read  (until  the  end-of-file 
is  reached).  Any  variables  on  the  new  PARAM  card  will 
have  the  new  values  assigned  to  them,  all  other  vari¬ 
ables  will  retain  their  previous  value  and  the  entire 
calculation  Is  repeated.  The  number  of  plots  which 
might  be  produced  by  such  a  run  is  equal  to  the  number 
of  PARAM  cards  in  the  input  deck. 

NOTE :  Only  one  set  of  data  points  for  Ncase®^  can 
be  used  within  a  given  run,  and  the  CASE4  card  and 
associated  NVAL,  YVAL  data  cards  need  only  be  supplied 
once  -  after  the  first  PARAM  card  which  contains 
Ncase=4 . 
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VII.  OUTPUT 
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1.  Printout 

a.  If  Ncase=4,  the  values  of  the  CASE4  variables  and  the 
associated  table  of  N(Y)/Nm  vs.  depth  values  is  printed. 

b.  For  each  parameter  set  the  values  of  all  variables  in 
the  PARAM  namelist  group  is  printed. 

c_.  For  each  parameter  set,  a  table  of  output  values  follows, 
giving  for  each  of  47  values  of  w: 

(1)  The  value  of 

(2)  The  value  of  the  integrand  INTI  ■  Sj_j  •  oi/CONSTl. 

(3)  NPT1  =  the  number  of  intervals  required  to  obtain 
the  integral  INTI. 

(4)  ACC1  -  the  relative  accuracy  achieved  in  the 
SIMPSON  Integration  INTI. 

(5)  The  value  of  Sji*.  The  column  is  headed  by  the  no¬ 
tation  REAL  or  IMAG  according  to  whether  the  values 
are  real  or  imaginary. 

NOTE :  If  ICALCC=0,  the  above  values  will  all  be 

zero . 

(6)  The  value  of  the  integrand  INT2  =  ^kJlmn  '  w/2C0NSTl. 

(7)  NPT2  =  number  of  intervals  required  to  obtain  the 
integral  INT2. 

(8)  ACC2  -  the  relative  accuracy  achieved  in  the  SIMPSON 
integration  INT2. 

(9)  The  value  of  Sicomn,} .  The  notation  REAL  or  IMAG  indi¬ 
cates  whether  the  values  are  real  or  imaginary. 

NOTE:  If  ICALCC=0,  the  four  preceding  values  will 
all  be  zero. 


The  values  of  Sij  and  Sk&mn  decrease  as  frequency  ( o)/2tt )  in¬ 
creases.  There  is  a  cutoff  value  VALMIN  such  that  once  the 
9  spectrum  being  calculated  is  less  than  VALMIN,  the  computation 

for  the  spectrum  ceases.  The  value  is  now  set  at  (10)-3,  but 
could  easily  be  changed. 
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The  main  reason  for  printing  the  NPT,  ACC  figures  was  to 
check  on  the  behaviour  of  the  function.  If  an  output  appeared 
with  NPT1  or  NPT2  ■  512,  it  would  Indicate  that  the  accuracy 
requested  had  not  been  reached  and  the  iteration  had  termi¬ 
nated  because  doubling  the  number  of  intervals  would  exceed 
the  limit  of  1000.  In  this  case,  further  examination  would  be 
indicated  to  see  if  some  error  had  been  introduced. 

2.  Plots 

When  the  plot  option  is  used,  each  run  will  produce  either 
one  or  two  plots  for  each  parameter  set  in  the  run  (S^  or 
Skimn  or  both).  The  Plot  is  1°"  by  10"  on  6-cycle  log  paper 
the  abscissa  being  frequency  (w/2iT«hz)  which  varies  from  0  to 
1.  The  ordinate  scale  is  chosen  in  such  a  way  that  the  maximum 
computed  value  of  the  function  (S^j  or  Sk^mn)  will  be  included 
on  the  plot,  but  in  no  oase  are  values  of  the  ordinate  less 
than  (10)-3  plotted.  Each  plot  is  labeled  completely  with 
the  values  of  all  the  parameters,  plus  an  alphanumeric  des¬ 
cription  of  the  VAIS ALA  frequency  assumption  used. 
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ADDENDUM  1 


TRIGONOMETRIC  POLYNOMIALS 


Let* 

X(w,(J>d)  “  cos2w  cos2<^d  +  sin2(j>i3 
Then  the  G.^  functions  are  defined  by: 

Gil (w )  "  ?  cos2(w-a)  X (w ,4>d> 

°22(W)  “  T  X  ( vr , ) 

g33(w)  -  ^  sin2  (w-a)  \(vr,<PD) 

°12  "  a2l^w^  "  IT  cos(w-a)  X(w,*D) 

Gig  ■=»  a3i(w)  =■  ^  cos(w-a)  sln(w-ot)  X(w»<{>D) 

G32(w)  -  G23  (w )  -  ^  sin(w-a)  X(w,<J>d) 
and  the  Gk^mn  functions  by: 

Gi2i2(w)  “■  £  cos2 (w-a )  X(wt<|>D)  =  \  Gn(w) 

Gi^is  (w ;  =  sin22  (w-a )  X  (w ,  4>d  ^  ^  sln2  2(w-a)  G22(w) 

G232o(w)  -  £  sin2(w-a)  X(w,<j>D)  B  |G33(w) 

Gl223(w)  =  ^  sin(w-a)  cos(w-a)  X(w,<()p)  “  ^  sin(w-a)  Gi2(w)/i 
G123l(w)  -  sin  2 (w-a)  cos(w-a)  X(w,<t>D)  -  sin  2(w-a)Gi2(w) 
G2313'w')  =  “T^  sln  2(w"a)  sin(w-a)  X  (w ,  sin  2(w-a)G32(w) 


The  Gij  used  here  are  denoted  by  the  lower  case  letter  in  [1], 
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ADDENDUM  2 


CASE  4,  SET  1  PARAMETERS  AND  DATA  POINTS 


The  STD  data  used  for  Ncase=4  runs  was  obtained  from  a 
ourve  of  N(Y)/Nm  v3.  Y  (See  Figure  1).  Fifty  points  were 
selected  from  the  curve  in  such  a  way  that  a  good  approxima¬ 
tion  to  the  curve  could  be  obtained  by  connecting  the  points 
with  straight  lines.  The  pairs  of  values  used  are  as  follows 


n 

N(Y?/Nm 

=1 

n 

NUIZNm, 

££ 

1 

.3 

0 

26 

.  880 

65.0 

2 

.325 

5.0 

27 

.882 

70.0 

3 

.350 

8.5 

28 

.880 

73.0 

4 

.375 

10.5 

29 

.675 

75.0 

5 

.40 

12.2 

30 

.870 

77.0 

6 

.5 

16.5 

31 

.855 

80.0 

7 

.  6 

19.0 

32 

.810 

85.O 

8 

.850 

26.5 

33 

.770 

90.0 

9 

.930 

30,0 

34 

.730 

94.0 

10 

.960 

31.5 

35 

.700 

98.0 

11 

.975 

32.5 

36 

.625 

110.0 

12 

.990 

34.0 

37 

.560 

120.0 

13 

.995 

34.5 

38 

.500 

130.0 

14 

1.0 

36.7 

39 

.440 

140.0 

15 

.990 

39.0 

4o 

.390 

150.0 

16 

.980 

40.0 

41 

.360 

156.0 

17 

.965 

42.0 

42 

.330 

164.0 

18 

.950 

44.0 

43 

.  300 

174.0 

19 

.930 

47.5 

44 

.280 

181.0 

20 

.915 

50.0 

45 

.260 

188.0 

21 

.910 

51.0 

46 

.250 

193.5 

22 

.900 

52.5 

47 

.240 

198.0 

23 

.890 

56.0 

48 

.225 

207.0 

24 

.885 

57.5 

49 

.220 

215.0 

25 

.880 

62.0 

50 

.210 

230.0 

The  parameters  associated  with  this  set  of  data  points 


are : 


Nset  =  1 
NPT  =>50 


Nm  -  .02182 
NMIN  -  .21 
YMIN  -  -230. 

Cp  -  .004756 
P  -  1.5 

so  that  for  depth  Y  <  -230,  N(Y)  -  (. 02182) (. 21)e ^Y+230^/130° 
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ADDENDUM  3 


NAMELIST  INPUT  CARDS* 

When  a  group  of  variables  is  identified  by  a  NAMELIST 
group  name,  input  and  output  of  these  variables  is  possible 
using  unformatted  data  cards  in  which  any  subset  of  the  vari¬ 
ables  in  the  group  may  be  assigned  values  as  follows. 

The  card  must  contain  a  $  in  Column  2  (Column  1  left  blank) 
immediately  followed  by  the  NAMELIST  group  name  with  no  em¬ 
bedded  blanks.  After  that  follows  a  string  of  data  items  of 
the  form  variable  *  constant,  (Example:  A-1.57,  kk“5»  CON  * 
22.3E-5)  separated  by  commas,  the  last  record  terminating  with 
a  $  instead  of  a  comma. 

More  than  one  card  can  be  used  to  oontain  a  NAMELIST  in¬ 
put  set  provided  that  the  first  column  of  each  card  is  left 
blank,  and  that  each  card,  except  the  first  one  ends  with  a  con¬ 
stant  followed  by  a  comma. 

To  illustrate,  if  the  following  set  of  NAMELIST  input 
cards  were  submitted  for  an  INTWAVE  run,  with  7  data  cards  con¬ 
taining  the  50  pairs  of  data  points  NVAL,  YVAL,  immediately 
after  the  CASE4  input  card,  the  run  would  produce  three  plots: 

1.  For  CASEit,  Set  1  data  with  a-  -  <fi,  V-100,  Y-50, 

p-1,5 

2.  For  CASE3  with  Ym=50,  a-  ir/2,  *d  -  0,  V-100,  Y-50, 

P“1 

3.  Same  as  (2.)  except  that  p-2. 


*This  description  applies  specifically  to  the  0006400  computer 
using  the  NOS/BE  operating  system. 
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NOTE:  On  any  PARAM  card  only  those  variables  whose  values  are 
to  be  changed  need  be  assigned  values.  The  variables  in  the 
group  which  are  not  assigned  values  will  be  taken  as  zero, 
provided  a  CLEAR  command  is  used  in  the  control  card  deck. 


ADDENDUM  4 


CORRELATION  COEFFICIENTS 


The  INTWAVf:!  program  can  also  be  made  to  generate  the 
following  three  sets  of  correlation  coefficients: 


f*1231<“) 

pl232(u)) 

pl323(ai) 
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where  the  quantities  are  defined  by  Equation  5,  all 

parameters  but  the  subscripts  being  kept  constant  over  the  set 
of  calculations. 

In  order  to  compute  these  correlation  coefficients  in  a. 
run  of  INTWAVE,  the  run  must  be  set  up  in  the  following  way: 

a.  A  new  option  parameter,  ICORR  must  be  added  to  the 
OPTIONS  namelist  input  card,  with  value  ■  1.  The 
value  of  ICALCG  must  also  equal  1. 

b.  The  first  PARAMS  namelist  Input  card  must  define  all 
the  required  parameters  including  one  of  the  six  sets 
of  kfcmn  subscripts  required  for  the  correlation  co¬ 
efficients. 

c .  The  next  five  PARAM  namelist  input  cards  must  change 

only  the  values  of  kk,  ll,  mm,  and  nn  in  any  order 
such  that  the  six  required  arrays  are  calculated. 
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d.  Correlations  for  other  parameter  combinations  can 
follow  the  first  set  of  calculations  by  repeating 
the  entire  sequence  of  six  PARAM  input  cards,  the 
first  one  in  each  set  defining  the  changed  para¬ 
meter  (s  )  . 

To  illustrate,  for  one  set  of  correlation  coefficients, 
assume  the  first  PARAM  card  contains  li"l,  jj»2,  kk»l,  JU»2. 
Then  a  set  of  five  PARAM  cards  which  would  produce  the  re¬ 
quired  arrays  might  be: 


1.  IPARAMS 

2.  $PARAMS 

3.  $PARAMS 

4.  $PARAMS 

5.  $PARAMS 


LLe3 ,  NN-3  $ 

KK*2 ,  MM" 2  $ 

MM“1$ 

KK-1,  LL"2,  MM"3 ,  NN-1  $ 
NN“2  $ 


Note  that  since  the  i3jc^mn  calculations  were  terminated 
when  the  values  became  less  than  VALMIN"10-^,  the  correlation 
will  be  given  as  0  if  any  of  the  spectrum  values  in  the 
numerator  <  10“3,  Also  if  any  of  the  speotrum  values  in  the 
denominator  <  10-3,  the  correlation  value  iB  undefined. 


After  the  six  required  skAmn  arrays  have  been  computed 
and  printed,  the  correlation  coefficients  will  be  computed 
and  printed  for  the  47  values  of  w  used  -  the  table  giving 
values  of  m/27r,  £1232^^'  and  p1323*'a)^' 
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